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Using a super-operator formulation of linearized time-dependent density-functional theory, the dynamical 
polarizability of a system of interacting electrons is given a matrix continued-fraction representation whose 
coefficients can be obtained from the non-symmetric block-Lanczos method. The resulting algorithm allows for 
the calculation of the full spectrum of a system with a computational workload which is only a few times larger 
than that needed for static polarizabilities within time-independent density-functional perturbation theory. The 
method is demonstrated with the calculation of the spectrum of benzene, and prospects for its application to the 
large-scale calculation of optical spectra are discussed. 
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The past two decades have witnessed the tremendous suc- 
cess of density-functional theory (DFT) f?, '5] in describing 
and predicting various properties of systems of interacting 
electrons, ranging from atoms and molecules to solids and 
liquids. In its original formulation, DFT only applies to the 
electronic ground state. This limitation was lifted by Runge 
and Gross (RG) |0] who generalized DFT to time-dependent 
systems. According to the RG theorem, for any given (t = 0) 
initial state of the many-body system, the external, gener- 
ally time-dependent, potential acting on it is uniquely de- 
termined by the time evolution of the one-electron density, 
n{r,t), for t > 0. Based on this theorem, it is possible to 
formally establish a time-dependent Kohn-Sham (KS) equa- 
tion from which various one-particle properties of the system 
can be obtained as functions of time. The resulting theoretical 
framework is usually referred to as time-dependent density- 
functional theory (TDDFT). Unfortunately, little is known 
about the exchange-correlation (XC) potential in ordinary 
DFT, and even less is known about it in the time-dependent 
case. Most of the existing applications of TDDFT are based 
on the adiabatic local density (ALDA) or adiabatic gener- 
alized gradient approximations, which amount to using the 
same functional dependence of the XC potential upon density 
as in the static case. In spite of the crudeness of these approx- 
imations, the optical spectra calculated using them is in some 
cases almost as accurate as those obtained from more compu- 
tationally demanding many-body approaches |4]. TDDFT is 
in principle an exact theory. Progress in understanding and 
characterizing the XC functional will thus substantially in- 
crease the predictive power of TDDFT while keeping its com- 
putational requirements at a significantly more modest level 
than methods based on many-body perturbation theory. 

Linearization of TDDFT with respect to the strength of 
some external perturbation to an otherwise time-independent 
system leads to a non-Hermitean eigenvalue problem whose 
eigenvalues are excitation energies of the system, and whose 



eigenvectors can be used to calculate appropriate oscilla- 
tor strengths f3\. Not surprisingly, this eigenvalue problem 
has the same structure as that arising from time-dependent 
Hartree-Fock theory j^, and the dimension of the result- 
ing matrix (the Liouvillian) is thus twice the product of the 
number of occupied states, N^, with the number of empty 
states. The diagonalization of such a large matrix can be 
accomplished using iterative techniques, often in conjunc- 
tion with the so-called Tamm-Dancoff approximation (TDA) 
which amounts to enforcing Hermiticity by neglecting the 
anti-Hermitean component of the Liouvillian fr]. Many of 
the existing molecular applications of TDDFT have been per- 
formed within such a framework which suffers however from 
two major limitations. The first is that many individual eigen- 
values/eigenvectors must be calculated in order to model the 
spectral properties in any given frequency range, and the num- 
ber of such eigenpairs increases with the size of the system. 
The second limitation is that a straightforward implementa- 
tion of this method preliminarily requires the full diagonaliza- 
tion of the unperturbed Kohn-Sham Hamiltonian, a task which 
may become prohibitive for large systems and/or when the 
one-electron basis set is very large (as may be the case with 
plane waves or real-space grids). The first of these problems 
is sometimes dealt with by directly calculating the relevant re- 
sponse function(s), rather than individual excitations The 
price paid in this case is the manipulation (inversion, mul- 
tiplication) of large matrices for any individual frequency, a 
task which may again be impossible for large systems/basis 
sets. Even so, the full diagonalization of the unperturbed KS 
Hamiltonian cannot be avoided. 

An alternative approach to TDDFT, first proposed by Ya- 
bana and Bertsch [8], avoids diagonalization altogether. In 
this approach, the TDDFT KS equations are solved in the time 
domain and susceptibilities are obtained by Fourier analyz- 
ing the response of the system to appropriate perturbations in 
the linear regime. This scheme has the same computational 
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complexity as standard time-independent iterative methods in 
DFT. For this reason, real-time methods have recently gained 
popularity in conjunction with the use of real-space grids |3, 
and a similar success should be expected using plane-wave 
(PW) basis sets. The main limitation of these methods is that 
the time step needed for a stable integration of the TDDFT KS 
equations is small (of the order of lO^'^ fs in typical pseudo- 
potential applications) and decreases as the number of PW's 
(or real-space grid points) increases. 

In this Letter we propose a novel way to calculate opti- 
cal spectra in the frequency domain — thus avoiding any ex- 
plicit integration of the TDDFT KS equations — which does 
not require any diagonalization (of either the unperturbed 
KS Hamiltonian, or the TDDFT Liouvillian), nor any time- 
consuming matrix operations. To this end, we first express a 
generalized susceptibility as a matrix element of the resolvent 
of the Liouvillian super-operator, defined in some appropriate 
operator space, between suitable operator states. This matrix 
element is then evaluated using a Lanczos continued-fraction 
technique. Each link of the Lanczos chain — that is calculated 
once for all frequencies — ^requires a number of floating-point 
operations which is only a few times larger than that needed 
by a single step of the iterative calculation of a static polariz- 
ability within time-independent density-functional perturba- 
tion theory (DFPT) floHll. 

Let us first consider an external perturbation whose Fourier 
transform will be indicated by X{lu)V' ext{^,^), where X{ljj) 
conventionally indicates the strength of the perturbation. Let 
(p° (r) and be the ground-state valence KS orbitals and en- 
ergies, and Lpy{r,t) the solution of the TDDFT KS equation 
with the initial condition (pv{j^, 0) = ¥'^(r). We define the or- 
bital response functions as (p'^{r,t) = e"''*(iy9„(r, i) — (y5°(r)), 
and indicate their Fourier transforms as (p^{r) = (^^(r,w), 
and (Py{r) — ffiyiv, —w). The ipf orbitals can be chosen to 
be orthogonal to the Kohn-Sham occupied manifold, and they 
can be easily shown to satisfy the coupled linear equations: 

Lo^tir) = (iJ^5-e„)^+(r)+AF'(r,c.)^°(r), 
~iu^-{r) = (H'ks~^v)^-{v)+PcV\v,u)^l{v),{\) 

where H'^c; is the ground-state KS Hamiltonian, Pc is the 
projector onto the KS empty-state manifold, V'{y,lo) — 
Vl^i{Y,Ld) + J K(r, r')n'(r', a;)(ir' is the perturbing KS po- 
tential, n'(r, = + 'Pvi'^)) is the lin- 
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ear variation of the electron density, and ^(r, r') = \T-r'\ 
^^Tpy- is the Hartree plus XC kernel. In Eq. Q 

^ ' n—n° 

use has been made of the fact that the perturbation is real 
[V'*{y, —uj) = V'{r,oj)], while in the expression for n', 
the (^°'s are assumed to be real due to time-reversal invari- 
ance. Following Refs. and 1 1 si we now consider the super- 
vector \X) = |x,y}, where X and y are batches of orbitals, 
X = {xy{r)}, y = {yv{r)}, defined as: x„(r) = 5[<y3+(r) + 
(p~ir)], and yy{r) = ^[(p+{r) - "^"(r)]. The fundamental 



TDDFT equations, Eq. Q, can then be formally written as: 

(c^-/:)|x,y) = |0,v), (2) 

where the batch V is defined as V = {PcV^^f{r, uj)(pl{r)}, the 
Liouvillian sMper-operafor £ is defined as: £|x,y) = \V ■ 
y, {V + K.) ■ x), and 2? and JC are Hermitean operators acting 
on batches of response orbitals, U = {uv{y)}, as: 

V-U = {(7?^s-e„K(r)}, (3) 
IC-U = I^SWE / ^^(r,r')¥':.(r'K'(r')dr'|. (4) 

Eq. (|2} gives the response of the system to the external pertur- 
bation, as a function of frequency. When V^^^ = 0, this equa- 
tion reduces to a (non-Hermitean) eigenvalue problem whose 
eigenvectors describe the free oscillations of the system cor- 
responding to electronic excitations. This is essentially equiv- 
alent to Casida's formulation of TDDFT |5]. 

In practice, one is seldom interested in the response of the 
system to the most general perturbation, or in individual exci- 
tation energies, but just in the frequency-dependent response 
of some specific property to some specific perturbation. Let 
us consider an observable. A, whose time-dependent linear 
response is given by: A{t) — 2Ke'^^{(f°\(p[,{t))- Assuming 
that A is time-reversal invariant, the Fourier transform of A{t) 
can be written as: A{w) = E,{{v°\A\ip+) + {ip°\A\ip-)) = 
2(a, 0|x,y), where a — {Aip",}. A generalized susceptibility, 
XAv{^), can be defined as the derivative of A{lu) with respect 
to the strength of the perturbation, X{u!). Using Eq. (|2j and the 
linearity of A{lu) with respect to X{u!), the susceptibility can 
be written as: 

XAvM = 2(a,0|(co-/:)-i|0,v}. (5) 

The results obtained so far and embodied in Eq. (|5j can 
be summarized by saying that within TDDFT any general- 
ized susceptibility can be expressed as an appropriate off- 
diagonal matrix element of the resolvent of the Liouvillian 
super-operator. In the following we show how such a matrix 
element can be conveniently calculated using a matrix contin- 
ued fraction resulting from the non-Hermitean block-Lanczos 
algorithm (NHBLA) 114). 

Let us define a block, |Q), as a pair of orthogonal super- 
vectors: |Q} = {IQi), IQ2}}- The scalar product between 
two blocks, = (P|Q} is defined as the 2x2 matrix: Sij — 
(PilQj), and the action of a super-operator on a block is de- 
fined as the block whose elements are the result of the ac- 
tion of the super-operator on each of the two elements of the 
original block: C\{Qi,Q2}) = \{CQi,CQ2})- Given two 
starting blocks, |Q^) and |P^} such that = Sij, the 

NHBLA allows one to generate a sequence of block pairs, 
{|Q"},|P")}, such that: {Pp\Qf) = SmrAj, and C = 
where: T™" = {PmQT) = 
afjSmn + bfjSm+i.n + c^(5m,n+i is a block-tridiagonal ma- 
trix. Using these relations, the projection of the resolvent of 
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FIG. 1 : Absorption spectrum of benzene calculated using the Lanc- 
zos method with different numbers of recursion steps: 1000 (plum), 
2000 (red), 3000 (green), and 6000 (black). The inset compares the 
6000-step spectrum (black) with that obtained using the real-time 
propagation method (orange). 



the Liouvillian over the starting block pair can then be easily 
expressed as a matrix continued fraction: 

(Pi|(a;-Z:)-i|Qi) = J , (6) 

w - ai + b2 C2 

uj — + ■ ■ ■ 

where the a's, b's, c's, are 2 x 2 matrices. If the starting block 
is chosen as \Q\) = |a, 0), and \Ql) ~ |0, v), the generalized 
susceptibility, Eq. Q, is then the (1,2) matrix element of the 
2x2 matrix given by Eq. (|6}. Without going into the details 
of the NHBLA, suffice it to say that its implementation does 
not require the explicit calculation of the Liouvillian super- 
operator, nor even of the unperturbed KS Hamiltonian, but 
just the availability of a black-box computer routine which, 
for any given batch of response functions, U ~ {wi,(r)}, re- 
turns 2?|u) and /C|u), according to Eqs. (|3} and Each step 
of the NHBLA essentially involves two calls to such a rou- 
tine whose computational cost is roughly the same as that of a 
single iteration in a static DFPT calculation. 

The theory described above was implemented on top of the 
PWscf plane-wave pseudopotential code which is part of the 
Quantum ESPRESSO distribution 1 15|. As a test to demon- 
strate our methodology, we have calculated the absorption 
spectrum of benzene, a system for which many excited-state 
calculations already exist [ [isl Holl . some of which were per- 
formed within TDDFT 11911 . The KS equations for an iso- 
lated benzene molecule were solved using periodic bound- 
ary conditions (PRC) in a tetragonal simulation cell whose 
base, parallel to the carbon ring ('a;?/' plane), has a side length 
30.0 A, and with height 2/3 of this; the CC and CH bond 
lengths were set to 1.40 A and 1.09 A, respectively. A sim- 
ple LDA exchange-correlation functional was adopted IT^ . 
whereas norm-conserving pseudopotentials from the PWscf 
table were used I17i] with a PW kinetic-energy cutoff of 60.0 



Ry. The PBC calculation of dipole matrix elements was per- 
formed as explained in Sec. C.2 of Ref. iljj. The absorption 
coefficient was obtained as: /(w) oc a;Imx(w), where x{^) 
is the electric dipole susceptibility, Eq. (|5}, calculated at com- 
plex frequencies with an imaginary part of 0.27 eV. 

The convergence properties of our block-Lanczos algo- 
rithm are displayed in Fig.^where we report the absorption 
spectra of benzene as calculated for light polarized in the xy 
plane, for different numbers of recursion steps. We see that 
2000-3000 steps are sufficient to ensure convergence for en- 
ergies up to « 15 cv. Our most converged spectrum (corre- 
sponding to 6000 recursion steps) is then compared with that 
obtained using the real-time propagation method (see the in- 
set). The agreement between the spectra obtained with these 
rather different methods is practically complete. Although a 
thorough theoretical analysis of the convergence properties of 
our algorithm is beyond the scope of the present letter, we 
would like to point out a few facts that, we believe, will de- 
serve further attention. We first notice that, not unexpect- 
edly, the convergence properties deteriorate with increasing 
frequency: the lower the frequency, the better the conver- 
gence. Second, the convergence rate of the algorithm seems 
to be somewhat affected by the condition number of the Liou- 
villian which, in turn, depends on the size of the PW basis set: 
the higher the PW kinetic-energy cutoff, the worse the con- 
vergence 1 20]. Furthermore, the numerical instabilities which 
are known to plague Lanczos diagonalization algorithms lll4ll 
seem to have little, if any, effect on the calculation of the re- 
solvent matrix elements through Eq. (|6|l. As far as we can say, 
Eq. (|6j can be pushed as much as needed to reach any desired 
level of accuracy. Finally, the non-Hermitean character of the 
Liouvillian seems to affect somewhat the efficiency of the al- 
gorithm, as can be seen from the performance of the TDA that 
we examine now. 

In Fig.Hwe compare the experimental absorption spectrum 
of Ref. 1 21 1 with those obtained from our TDDFT method, 
with and without use of the TDA. The agreement between the 
calculated and experimental spectra is very good, as was al- 
ready known from previous TDDFT calculations for benzene 
lll9ll . The quality of this agreement is due to the nature of 
the molecular orbitals involved in the transitions which domi- 
nate the low-lying part of the spectrum (tt, tt*, and, to a lesser 
extent, cr), which are little affected by the wrong asymptotic 
behavior of the ALDA XC potential. The absorption spectra 
reported in Fig.|2j though resulting from an average over dif- 
ferent polarizations, are dominated by the (tt* ^ tt) 
transition which is only allowed when the light is polarized in 
the plane of the molecule. This transition is mainly responsi- 
ble for the first strong absorption peak experimentally found at 
6.94 eV, and predicted by TDDFT at 6.83 eV. It is interesting 
to notice that the z component of the spectrum displays a weak 
peak at 6.55 eV, which is not visible in Fig.|2](its intensity is 
more than 10 times smaller than the xy peak), and which cor- 
responds to a (tt* ^ cr) ^A2u transition. In the independent- 
electron approximation, this transition would have a higher 
energy than which corresponds to the HOMO-LUMO 
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FIG. 2: Comparison between the converged Lanczos spectrum ob- 
tained with the full non-Hermitean Liouvillian (black), the Tamm- 
Dancoff approximation (red) and the experimental results (light 
blue). The inset shows the convergence of the TDA spectrum with 
respect to the number of recursion steps: 250 (blue), 500 1000 (ma- 
genta), and 2000 (orange). 



gap. The red shift of the ^A2u transition is therefore due to 
the effects of the electron-electron interaction which are ap- 
proximately accounted for in ALDA-TDDFT. The ^A2u tran- 
sition has never been detected directly in absorption experi- 
ments, but its existence (as well as its location near, possibly 
at a lower frequency than, the strong ^Eiu transition) was in- 
ferred from Raman scattering experiments 1E2I1 . The proxim- 
ity of the ^Eiu and ^A2u, as well as the much smaller intensity 
of the latter, were also confirmed by accurate coupled-cluster 
calculations As far as we can tell, it is not impossible 

that the little shoulder observed at 6.19 eV in the absorption 
spectrum |21] and attributed to a vibron-assisted forbid- 
den excitation, is actually due to a weak ^A2u allowed transi- 
tion. 

Use of the TDA does not change much the overall ap- 
pearance of the spectrum, nor the positions of the peaks, the 
main difference regarding their intensities. It is worth notic- 
ing that the convergence of TDA calculations appears to be 
much faster than when using the full non-Hermitean form of 
the Liouvillian (see the inset of Fig.|2}. 

We believe that the method presented in this letter, while 
not touching our still serious ignorance about the form of 
time-dependent XC functionals, will open the way to a sys- 
tematic study of systems which are too large to be treated with 
currently available methods. The already rather favorable nu- 
merical features of our method will be further improved either 
by implementing the use of ultra-soft pseudopotentials (which 
will allow reduction of both the size of one-electron basis sets 
and the condition number of the Liouvillian), and by devis- 
ing optimal strategies for restarting the Lanczos chain using 
approximate schemes. Work is in progress along these lines. 
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